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I. INTRODUCTION 



Modeling temporal fluctuations in avalanching systems 
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We demonstrate how to model the topphng activity in avalanching systems by stochastic differ- 
ential equations (SDEs). The theory is developed as a generalization of the classical mean field 
approach to sandpile dynamics by formulating it as a generalization of Itoh's SDE. This equation 
contains a fractional Gaussian noise term representing the branching of an avalanche into small 
active clusters, and a drift term reflecting the tendency for small avalanches to grow and large 
avalanches to be constricted by the finite system size. If one defines avalanching to take place when 
QQ , the toppling activity exceeds a certain threshold the stochastic model allows us to compute the 

f^ ■ avalanche exponents in the continum limit as functions of the Hurst exponent of the noise. The 

t"*^ ' results are found to agree well with numerical simulations in the Bak-Tang-Wiesenfeld and Zhang 

^\] , sandpile models. The stochastic model also provides a method for computing the probability density 

functions of the fiuctuations in the toppling activity itself. We show that the sandpiles do not belong 
to the class of phenomena giving rise to universal non-Gaussian probability density functions for 
the global activity. Moreover, we demonstrate essential differences between the fluctuations of total 
fvj . kinetic energy in a two-dimensional turbulence simulation and the toppling activity in sandpiles. 

(N 

o 

■^^ ' The aim of this paper is to present a consistent framework for the modelling of temporal fluctuations, including 

^ ' definition and computation of avalanche exponents, in sandpile (height) models such as the Bak-Tang-Wiesenfeld 

• ^ (BTW) and Zhang models [l|, [^ . One of the defining properties of self-organized criticality is that avalanche duration 

and size are quantities subject to scaling, i.e. pdm{T) ~ t~" andpsizc(s) ~ s~^ I2l,lj]- The calculation of the exponents 

a and v in the thermodynamic limit N —^ oo {N'^ is the number of sites in the d-dimensional lattice) has proven to 

be a difficult task in dimensions d — 2 and d = 3. This is partly due to the lack of simple finite-size scaling in models 

such as the BTW model [3]. Moreover, it has recently been pointed out [5'] that the difficulty may originate from 

\^^ ' the fact that the r and s do not scale when defined in the traditional sense, which is to consider the duration of an 

7-H , avalanche as the time interval where toppling takes place between two successive zeroes in the toppling activity. In 

^^ this paper we shall denote avalanches defined this way as type-I avalanches. 

CO The scaling property can be restored, however, if one defines the duration of an avalanche as the time interval 

J — . when the toppling activity x{t) exceeds a prescribed threshold Xth- We shall use the term type-II avalanches when 
(T^ ' the start and end of an avalanche is determined via such a threshold criterion. The idea of using a threshold on the 
OO , toppling activity in the definition of avalanches was first introduced in [6|, where it was argued that any avalanche 
^^ ■ analysis of real-world activity time series must define avalanches from a threshold, since there is no way to uniquely 
determine whether a non-negative continuous- valued experimental quantity is actually zero, or just small. For type-II 
avalanches it can be shown by numerical simulation that the quiet times between avalanches in the BTW model are 
. _ power-law distributed. For type-I avalanches the quiet times only depend on the statistics of the driver, which is 
J3 ■ usually assumed to be Poisson distributed. 

The present work represents the first systematic investigation of type-II avalanche statistics for the BTW and Zhang 
sandpiles. We are particularly interested in the continuum limit where the system size L = 1 is considered fixed and 
the spatial resolution increases as iV ^ oo. The power-law statistics of avalanche observables have cutoffs for large 
avalanches due to the finiteness of the system, but as we increase the resolution we see an increasing range of scaling 
for smaller avalanches. As A^ ^ oo we keep the threshold fixed relative to the time-averaged activity (x), which scales 
as (x) ~ iV^i, where < Di < 2. For the BTW sandpile, numerical simulations yields Di « 0.86 [^]. This means 
that the number of overcritical sites corresponding to the threshold value diverges like Xc ~ N^^ in the limit N —> oo. 
We model the toppling activity in the continum limit by a stochastic differential equation [7[ for a normalized 
toppling activity X(t). In its simplest form this equation is on the form 

dX{t) =(7^/X{t)dW{t), (1) 



where W{t) is the Wiener process. Without the factor y^X{t) on the right hand side we would simply have that 
X{t) — aW{t) + Xq is a Brownian motion with diffusion coefficient D = a^/2. This factor, however, gives rise to 
a non-uniform (X-dependent) diffusion coefficient D = a'^X{t)/2, and the stochastic process X{t) will have non- 
stationary increments. This model can be perceived as a continuous version of the classical mean field theory of 



sandpiles [g. We use its corresponding Fokker-Planck equation to derive that a = 2 and v = 3/2 when avalanches 
are defined in the type-I sense. This is the same results obtained by mean field theory [3, [j, Q. For type-II avalanches 
the effect of the non-uniform diffusion coefficient vanishes for avalanches of durations short compared to that of a 
system-size avalanche, and for the purposes of calculating a and v we can assume that the the toppling activity is 
a standard Brownian motion. Solving the Fokker-Planck equation for Brownian motion (or equivalently using the 
known distribution of first return times in Brownian motion) we obtain a = 3/2 and v = 4/3. 

For the modeling of non-trivial sandpile models (BTW and Zhang) the stochastic differential equation takes the 
form 



dX{t)^f{X)dt + ay/X{t)dWHit). (2) 

Two important generalizations of the stochastic model are included here. First, from numerical simulations of sandpiles 
we find that for small activities X{t) there is an effective positive drift term. This term is dominant for very small 
activity, since the diffusion term is negligible for very small X(t) due to the X-factor in the diffusion coefhcient. 
The positive drift term is X-dependent and quickly decreases as X increases, but strongly influences the avalanche 
statistics because it contributes to prevent avalanches from terminating when X{t) approaches zero. We believe that 
this effect is responsible for destroying the scaling of avalanche duration and size (when these are defined in the type-I 
sense). However, if the drift term is small compared to the diffusion term for X > Xth, the drift term will not affect 
the avalanche statistics if one employs a threshold Xth to define type-II avalanches. This explains why scaling of size 
and duration is restored when when type-II avalanches are introduced. Another generalization, which is essential for 
the avalanche statistics, is the Hurst exponent H of the noise term. The mean-field approach to sandpiles implicitly 
assumes that H = 1/2. However, this is not the case for the BTW and Zhang models. Actually, analysis of numerical 
simulations of the sandpiles show that H = 0.37 for the BTW model and H = 0.75 for the Zhang model. 

As in the mean-field model, the effect of the non- uniform diffusion coefficient vanishes as the threshold increases, and 
hence keeping the threshold Xth fixed and increasing N we can, for the purposes of computing a and v for avalanches 
where X never grows much greater than Xth, consider the toppling activity as a fractional Browian motion. Using the 
result of Ding and Yang [9|], that the first return time in fractional Brownian motion scales like ~ r^^^ , we obtain 
the general results a = 2 — H and v = 2/(1 -|- H). 

Although the drift term and the non-uniformity of the diffusion coefficient are not important to calculate the 
avalanche exponents for type-II avalanches whose duration are short enough not to be limited by the finite system 
size, they are important on the time scales where the toppling activity is a stationary process. These are scales 
sufficiently long that the toppling activity of avalanches is limited by the boundaries. The stochastic equation ([2|) 
is fully equipped to handle these time scales. A good example of the applicability of these aspects of the stochastic 
model is the computation of the probability density function (PDF) of the temporal fluctuations in the activity signal 
itself. 

For weakly driven sandpiles the PDFs of the fluctuations in the toppling activity are stretched exponentials. This 
result is reproduced by simulation of the stochastic model [Sj . For stronger driving the activity exhibits fluctuations 
which are more confined around a mean value where the drive and dissipation balance each other. It has been 
claimed that the PDFs of the toppling activity in sandpiles are examples of universal Bramwell-Holdsworth-Pinton 
(BHP) distributions [lO, [ll|) a certain class of asymmetric PDFs commonly seen in complex systems. Our sandpile 
simulations show that the BHP distributions can only be seen if one fine tunes the driving rate to a certain value, and 
for other driving rates the PDFs belong to a much wider class of distributions. For sufficiently strong drive Gaussian 
PDFs are observed. These can also be obtained from the stochastic model if one correctly models the drift term in 
this parameter range. 

The rest of the paper is structured as follows: In Sec. |TT] we explain and derive the stochastic model for the toppling 
activity. In Sec. IIIII we compute the avalanche exponents for type-I and type-II avalanches for the mean field case 
{H = 1/2) by solving a Fokker-Planck equation, and for H ^ 1/2 we compute the avalanche exponents for type- 
II avalanches as a function of H . The results allow us to predict the avalanche exponents for sandpile models by 
computation of the Hurst exponents. These results are then tested against numerical simulations of the BTW and 
Zhang models and are shown to agree well. In Sec. IIVI we present results which indicate that type-II avalanches 
exhibit so-called finite-size scaling, even though type-I avalanches do not. 

In Sec. |V] we use the stochastic theory to calculate the PDFs of the toppling activity signal, both for strongly 
driven sandpiles and in the weak driving limit. The method is finally applied to the fiuctuations in kinetic energy in 
a two-dimensional turbulence simulation. In this case the process is given by a different kind of stochastic differential 
equation: 

dX{t)=bdt + ce''^dW{t). (3) 

This equation gives rise to a Fischer-Tippet-Gumbel (FTG) distribution [13, [3 1 which is very close to the BHP. The 
differences between the stochastic models for sandpiles activity and kinetic energy fiuctuations in 2D turbulence may 
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FIG. 1: a) A realization of the toppling activity x{k) in the BTW sandpile. b) The increments Ax{k) — x{k + 1) — x{k) of the 
trace in (a), showing that Ax{k) is large when x{k) is large, c) Conditional PDFs oi x + Ax for x = 10, 20, 30 respectively, d) 
The conditional mean and variance of Ax versus x. 



represent an essential distinguishing feature between 2D turbulent dynamics and the kind of avalanching dynamics 
which are observed in the classical sandpile models. 
In Sec. IVIIwe summarize and conclude the work. 



II. THE STOCHASTIC MODEL 



Let x{k) denote the number of overcritical sites at time step A; in a sandpile. A common feature of (height-type) 
sandpile models is that the typical size of increments Ax is proportional to the square root of the toppling activity x. 
To be more precise, the conditional probability of an increment Ax(fc) = x{k + 1) — x{k), given x = x(k), is 



P{Ax\x) 



1 



\/2^ 



(4) 



In FiglT] this property is verified for the BTW model ( it holds in the Zhang model as well). This result can be 
explained as follows: At a given time k there are x = x(k) overcritical sites, which we can enumerate i = 1, 2, . . . , x. 
In the next time step fc — > A: + 1 each site i distributes energy to its neighbors, and will usually (always in the BTW 
model) become subcritical. If none of the neighbors receive sufficient energy to become overcritical, the contribution 
to Ax from site i is ^i = — 1. If exactly one of the neighbors become overcritical, then ^i = 0, and so on. For the 
two-dimensional models the maximal value of S,i is 4 (3 for the BTW model) since a site maximally can excite 4 
neighboring sites. Hence, for this case, we consider S,i to be random variables with realizations in {—1,0, 1,2,3,4}. 
The randomness originates from the local configuration in the vicinity of the overcritical site, which for these purposes 
is considered to be random. In other words, we think of the configuration on the lattice as a random background. 

As an approximation we consider the different realizations of S,i (at a fixed time k) as independent of each other. 
We also consider the distribution of S,i to be identical for all overcritical sites and for all times. In this approximation 
the total increment Ax can be written as a sum of independent, identically distributed, random variables Aa; = 
^1 + . . . + ^2, and by the central limit theorem we have @ provided that the local means (^) are zero. Then 
a^ = {^^) = (— l)^p_i + ... + 4^p4, where p = (p_i,...,p4) is the probability vector for the local increment 
processes. 

If the local processes at different times k are independent of each other, then x{k) is a Markov process which satisfies 
a stochastic difference equation 



Ax{k) = a^x{k) w{k) , 

where w{k) is a stationary, normalized, and uncorrelated Gaussian process, i.e. w{k) — W{k- 
is the Wiener process. Under a rescaling of time t = k6t and X{t) — x{k) 5x we have 



1)-W{k), where Wit) 



AX{t) ^ X{t + St)-X{t) = Sx(x{k + l)-x{k)) 



SxayG(k)w{k) ^ Sx^/"^ a ^/x{t) (w{k + 1) ~ W{k) 



/fe\i/2 



^) ay'X{f){w{t + 6t)-W{t)). 



\St 

In the last step we used the self-afRnity of the Wiener process, W{t/St)=St~^^^W{t). For 6x = St we have a well 
defined model in the limit St,5x — > 0, namely the Ito stochastic differential equation Eq. ([T]). 

The first generalization of this model is obtained if we relax the requirement that the local increment processes 
^i(fc) at time k are independent of the local increment processes ^j{k'), j = 1, . . . ,x{k') at previous times k' < k. In 
this case we need to model memory effects in the stationary Gaussian process 

Ax{k) 
w[k) — 



(Jy/x{k) 

From the power spectrum or the variogram of the activity signal from numerical simulation of the BTW and Zhang 
models we find that w{k) can be accurately modelled as a colored noise characterized by a Hurst exponent H . That 
is, w{k) — Wnik + 1) — Wnik), with Wh being a normalized (diffusion coefficient — 1) fractional Brownian motion 
(fBm). If we perform the rescaling t = kSt and X{t) = x{k) Sx in the case H ^ 1/2 we obtain 

^X{t) = [^y^\^X{r){wH{t + 5t)-WH(t)), 
and by requiring that 5x = 5t^ , with h = 2_ff, we obtain the stochastic differential equation 



dX{t) = (JyJX{t)dWH{t) . (5) 

If we assume that the the stochastic process X{t) is self-affine with self-affinity exponent ft-, i.e. X{st) = s^X{t), it 
is easy to verify that Eq. (O is invariant with respect to the transformation t —> st if ft, = 2H. Thus, the exponent 
h = 2H is the self-affinity exponent of the process X{t), where H is the Hurst exponent determining the color of the 
noise process driving the stochastic differential equation. Observe that the reason why h ^ H is the non-stationarity 
of the increment process due to the the factor ^/X(t) in Eq. ([5]). Note also that the case H — 1/2 corresponds to 
h = l. 

The self-affinity of X{t) described by Eq. ([5]) implies that there is no upper bound on the fluctuations on increasing 
time scales, i.e. Eq. ((S]) describes the activity of an infinite sandpile where the activity is never influenced by the 
system boundaries. From a physical viewpoint, however, it is more interesting to consider the dynamics of a flnite 
sandpile in the continuum (thermodynamic) limit, and for this purpose it is natural to let the scaling factor Sx depend 
on N such that Xn = xjfSx{N) is bounded in the limit N — > oo, for instance such that limjv_»oo max(X7v) = 1. 
Such a bound on X{t) can be obtained by the introduction of a drift term f{X) dt leaving the stochastic equation in 
the form of Eq. ^, where f{X) is negative for large X. The form of f{X) can be found from sandpile simulations 
by computing the conditional mean E{Ax\x) of the increments, and is shown in FigUJl. It appears that f{X) is a 
decreasing function, positive for small X and negative for large X, and f{X) = for a characteristic activity Xc ~ 1. 
Without the stochastic term the drift term establishes Xc as a stable fixed point for the dynamics. 

Since f{X = 0) > solutions of Eq. ^ with initial condition X{Q) > exist and are positive for alH > 0. This 
means that while Eq. ([5]) has solutions for which X{t) = after a finite time (avalanches terminate), avalanches 
described by Eq. ^ will never terminate in the meaning X{t) = (type-I avalanches). This signifies that type-I 
termination never occurs in the the continuum limit. If sandpiles of increasing N are simulated, and the type- 
I avalanche durations are computed in the rescaled coordinates, the durations generally grow without bounds for 
increasing N. The reason is that the effective threshold for type-I termination in a discrete sandpile is xn = 1, but 
in rescaled coordinates this threshold Xn — xn dx{N) goes to zero as A^ -^ oo. As this rescaled threshold vanishes 
the duration in rescaled time goes to infinity. Since the type-I termination is a discreteness effect, the resulting PDFs 
of avalanche durations depends on N (system discreteness) and are not power-laws. As we shall demonstrate later, 
the introduction of activity thresholds which are defined in the rescaled coordinates, and hence remain finite in the 
continuum limit, will give rise to PDFs of durations of type-II avalanches which converge to a specific power-law in 
this limit. 

Eq. ^ remains valid also for sandpiles which are driven by continuous feeding of sand during avalanches. The drive 
adds a positive contribution to the drift function f{X) for X < Xc, but a negative contribution for X > Xc, because 
for large activites f{X) mainly accounts for the increased boundary losses. The result is a steeper f{X), which tends 
to confine the activity closer to the fixed point Xc- On the other hand the increased drive also increases the diffusion 
coefficient (by increasing a) due to a larger number of new active clusters initiated per unit time. The net effect is 



a positive shift of X^ and that the fluctuations in X are confined to a smaller region around Xc- For sufficiently 
strong drive the range of variation in X{t) becomes so small that the diffusion coefHcient does not vary much, ft is 
nevertheless important to model it correctly in order to calculate the PDFs of the fluctuations in toppling activity. 

III. CALCULATION OF AVALANCHE EXPONENTS 

We denote the stochastic model Eq. ([T]) (which is Eq. ^ with H = f /2 and f{X) ~ 0) the mean field model of 
sandpiles. This is because its underlying assumptions and the results derived from it coincide with what is known 
as the mean field solution of sandpiles in the literature [8]. Eq. (jT]), together with its corresponding Fokker-Plack 
formulation can be used to calculate the avalanche exponents a and ly. The idea is that each avalanche corresponds 
to a realization X{t) with some initial condition X{0) = Xq {Xq ^ Xmax)- The avalanche propagates until the 
realization X(t) terminates a.t t — ti in the meaning that X{t) > ioi t < ti and X{ti) = 0. Calculating the 
ratio of surviving realizations at different times t in an ensemble will provide information about the distribution of 
avalanche durations. The avalanche size distribution can then be obtained by using a general relationship between 
the self-affinity exponent h = 2H = 1 and the duration statistics. 

The scenario outlined above can be mathematically formulated as follows: Let P{X, t) = Pr[X(t) exists and X{t) = 
X]. Then the the probability p{t) that an avalanche still runs after a time t (we call it the survival function) is given 

by 



/ 



pit)^ / PiX,t)dX, (6) 

Jo 

and the probability density function for durations is Pdur(T) = —p'{t). The density P{X,t) can be calculated by 
solving the Fokker-Planck equation 

dP a^ d^ 

on X G [0,oo), t € [0,oo), subject to an absorbing boundary condition limx^o ^^(^j ^) = ^^nd an initial condition 
PiX,0)=6iX-Xo). 

To correctly incorporate the absorbing boundary condition we let U = XP, and solve the corresponding Fokker- 
Planck equation for U with boundary condition U{0) = to get 

/•OO 

P{X,t)= G{X,Y,t)P{Y,0)dY, 

Jo 

where 

G{X,Y,t)^-sJ- / Ji(s\/F)ji(s\/X)exp( —tUds. 

Remark L The most elegant way to obtain the solution is to use the integral transform pair 

F{s) = - / F{X)Ji{sVx)VxdX 
2 Jo 

1 f°° 

F{X) = —= / F{s)Jl{s^/X)sds. 
yX Jo 

Taking the transform of the Fokker-Planck equation we obtain the ODEs 

— (S,i)= —P{s,t), 

which we can solve and take the inverse transform. It is also easy to solve the Fokker-Plack equation for U — XP by 
separation of variables. 



If P{X, 0) = S{X — Xo) the solution of the absorbing boundary problem is P(X, t) = G{X, Xq, t). Moreover 
,„ P,X.,) . ^ r,'MsV%) exp (-^.) .. ^ if e.p (-Mi 
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FIG. 2: a) Determining /i = 2H in the Zhang model through the calculation of the variance of the activity X{k) in avalanches 
still running at time t. On shorter time scales there is a Hurst exponent H = 0.5 (the first dashed line has slope 2h — AH — 2) 
on longer time scales the Hurst exponent is H = 0.75 (the second dashed hue has slope 2h — AH — 3). b) Determining h — 2H 
in the BTW model by the same method as in (a). In this case we have H = 0.37. 



From Eqs. 



(Ill), (111), and the absorbing boundary condition at X = wc find that 

2 Xq I 2 Xq 



dp 



— lim PiX, t) 

2 x^o ^ ' ' 



uH' 



exp 



aH 



Hence the PDF for avalanche durations r is 



Pdui(T) = -2"^ exp 



2Xo 



and for r >> 2Xo/cr^ we have pdm:[T) ^ t~", with a = 2. 

This result crucially depends on the correct formulation of the Fokker-Planck equation. If the stochastic process 
X{t) were a classical Brownan motion, the Fokker-Planck equation would have the form of the standard heat equation, 
and by performing the analogous calculations for this equation we get Pdui(''') ^ t~^/^. The same scaling relation 
('^ T^/^) is obtained if one (incorrectly) employs the Stratonovich formulation of the Fokker-Planck equation rather 
than the Itoh form [7] . 

From the derivation of the stochastic model we have seen that the process X{t) has a self- affinity exponent h = 2H 



(for H — 1/2 we have h — 1). This means that X{t) disperses like ' 
of the avalanche scales as 



t . For long avalanches this implies that the size 



X{t)dt 



t'^dt 



rh + l 



(9) 



This property is easy to check directly by studying the relation between duration and sizes of avalanches in sandpiles. 
We find that the relation holds for large avalanches, whereas for very small avalanches s ^ t^. If the survival function 
scales as p{t) ~ t^^ = t^"+^, then using ([9]) together with Psizo(s) ds = Pdui{T) dr yields that if Psize(s) ^ s^" , then 



l + h + 5 
1 + h 



1 



(10) 



In the case H — 1/2 {h — 1) and a — 2 this gives ly = 3/2. 

This mean- field solution is known to be quite correct for the random neighbor sandpile model [15J, but numerical 
simulation shows that is fails for type-I as well as type-II avalanches in the BTW and Zhang models. The computation 
of these exponents are shown for these models in Figures [2l [3l and [4] The computation of h presented in Fig. [2] for 
type-I avalanches yield H = h/2 = 0.37 for the BTW model, and H = 0.75 for the Zhang model. This invalidates the 
Fokker-Planck formulation, which can be strictly justified only for a white noise source term {H = 1/2). This is one 
reason for the failure of the mean-field approach, but there are also others. 

For avalanche duration and size type-I avalanches do not yield good power-law PDFs. The reason for this was 
discussed in the previous section: letting avalanches terminate when X]y{t) = in a sandpile with N"^ sites corresponds 
to using an effective threshold for termination in the rescaled coordinates X]^(t) which goes to zero as A^ ^- oo. The 
termination process depends on the discreteness of the system and one cannot expect convergence to scale-invariant 
behavior in the continuum limit. 
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FIG. 3: a) The survival function of avalanches in the BTW model computed with thresholds {X)/3 (type-II) and without 
thresholds (type-I, the dotted curves). The dashed line corresponds to a — 2 — H — 1.63. b) The probability density function 
of avalanches with size > s in the BTW model computed with thresholds {X)/3 and without thresholds (the dotted curves). 
The dashed line corresponds to i/ = 2/(1 + H) — 1.49. 
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FIG. 4: a) The probability of having avalanches of duration > r in the Zhang model computed with thresholds {X)/3 and 
without thresholds (the dotted curves). The first dashed line corresponds to a = 2 — H = 1.50 obtained with H = 0.50, and 
the second dashed line corresponds to a — 2 — H — 1.25 obtained with H = 0.75. b) The probability of having avalanches with 
size > s in the Zhang model computed with thresholds {X)/3 and without thresholds (the dotted curves). The first dashed line 
corresponds to i^ = 2/(1 + H) — 1.33 obtained with H = 0.50, and the second dashed line corresponds to // = 2/(1 + H) — 1.14 
obtained with H = 0.75. 



For type-II avalanches Figs. [3] and H] show good scaling for duration and size in BTW and Zhang models, but the 
scaling exponents differ from those of the mean-field approach. The discrepancy is partly due to the fact that the 
sandpile models have H ^ 1/2, but is also related to the observation that the introduction of a finite termination 
threshold Xj/j modifies the exponent h as it appears in Eqs. (0 and PH)) . In the following we shall demonstrate that 
it is possible to obtain analytical results for type-II avalanches from the stochastic model which are in agreement with 
the corresponding sandpile simulations. 

First we observe that omission of the drift term will only have effect on avalanches which are so large that they are 
strongly limited by boundary dissipation. Next, we notice that the effect of the positive drift term for small activities 
is eliminated for type-II avalanches if Xth > Xc- In other words, it is only the cut-offs of the power-law PDFs due to 
finite system size which are lost by this omission. Thus, by considering a threshold Xth which is not much smaller 
than the mean activity {X), and considering only avalanches which arc sufficiently small not to be influenced by the 
boundary dissipation, we have X(t) ~ Xth and thus can justify the substitution X[t) —f Xthit) on the right hand side 
in Eq. ^. This equation then reduces to the equation for an fBm with Hurst exponent H. For an_fBm the duration 



of avalanches is given by the return time statistics for Wh (t) , which is known to scale like t^ 



i.e. we have 



a = 2-H . 



(11) 



Now we need to address a slightly subtle point: the exponent h, as defined in Eq. ([9]), is h — 2H for times r so 
large that X{t) 3> -'^(O). Only on such time scales will the effect of the factor ^/xJt) in the stochastic term show 
up in the scaling. However, this makes sense only for type-I avalanches where we can choose X{0) <C {X). For 
type-II avalanches we have X{0) « Xth, and it is more natural to consider the opposite limit where r is so small that 
X{t) = X{t) — Xth ^ ^{0) ~ ^th- On these time scales the activity measured relative to the threshold level scales 
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like an fBm with Hurst exponent H, i.e. X{t) ^ t^ , because X{t) « Xth- The avalanche size of type-II avalanches 
defined as s{t) = L X{t) dt then scale as ^ r^+^. Thus, for type-II avalanches we have h = H, and hence Eq. (flOl) 
for this case reduces to 



1 + H 



(12) 



In the mean field limit H = 1/2 these exponents reduce to those given in the right hand column in Table HI These 
results obtained analytically by approximating the coefficient on right hand side in the stochastic equation by its 
threshold value has been verified by numerical solutions of the full equation. 

TABLE L Exponents in mean-field solutions [H = 1/2) of the stochastic equation with zero threshold (type-I) and large 
threshold (type-II). 





type-I avalanches 


type-II avalanches) 


H 


1/2 


1/2 


h 


1 


1/2 


a 


2 


3/2 


V 


3/2 


5/4 



Eqs. (fTTjl and p2|l show that the calculations of a and v reduces to determining the Hurst exponent of the normalized 
increment process 



w{k) 



Ax{k) 



which can easily be constructed from the toppling activity signal. The corresponding Hurst exponent of the motion 



W{k)=^w{k) 



4=0 

can be calculated using standard techniques such as taking the power spectrum or constructing variograms. Alter- 
natively we can find H by computing the variance of X{t) in an ensemble of realizations starting with small initial 
values of X. This variance scales like {X'^{t)) ~ i^'' = t^^ . Fig[2] shows this computation for the BTW and Zhang 
models. Since we find H = 0.37 for the BTW model the type-II scaling exponents are a = 1.63 and v = 1.49. This is 
verified by numerical calculation of a and v using a threshold (a;)/3, as shown in Figl3l 

In the Zhang model the process w(fc) is more complicated. Figl^] shows that W{k) has a Hurst exponent H = 0.5 
on short time scales and a different Hurst exponent H = 0.75 on longer timescales. This means that short avalanches 
should satisfy the mean-field solution a = 1.5 and i' = 1.33, whereas longer avalanches should have exponents a = 1.25 
and i^ = 1.14. All of these predictions are verified by the direct calculation of a and v as shown in FiglSl 

For increasing N the long-avalanche scaling dominates an increasing portion of the graph, so in the continuum limit 
the mean-field solution only prevails at infinitely small scales in rescaled coordinates. It is however a nice verification 
of our method to see that the relation between the avalanche exponents and the Hurst exponent correctly predicts 
the avalanche exponents on short time scales as well. 

These results on the BTW and Zhang models are summarized in table [Hi 

Remark 2. The numerical simulations of the Zhang model are run using the standard toppling rule z^ — > and 
Zj -^ Zi + Zi/A if Zi is overcritical and j is a nearest neighbor of i. Whenever the configuration has no overcritical 
sites a random site i is chosen with respect to uniform, probability and a mass e is added to this site: zi = Zi + e. In 
the simulations presented in this paper we use e = 0.1. In the strongly driven Zhang models presented in Sec. [B the 
feeding times (which can now be during avalanches) are Poisson distributed and we have used e = 0.25. 

For the BTW model we have used the standard toppling rule Zi —^ Zi — 4 and Zj = Zj + 1 through out the paper. As 
usual a mass e = 1 is added to a random site whenever the configuration is stable. 



TABLE II: Exponents for type-II avalanches computed from H = 0.37 in the BTW model, and H = 0.50 (short time scales) 
and H = 0.75 (long time scales) in the Zhang model. 





BTW 


Zhang 
(short times) 


Zhang 
(long times) 


H 


0.37 


0.50 


0.75 


a. 


1.63 


1.50 


1.25 


V 


1.49 


1.33 


1.14 
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FIG. 5: a) The moments (s'') plotted as functions of A*' for the size distribution of the BTW model ith respect to a threshold 



{X)/3. b) The shape of the structure function C(9)- The slope of the line is D 
is i^ — 1 = 0.49. This value is indicated by the dotted vertical line. 



2.81 and the intersection with the first axis 



IV. FINITE-SIZE SCALING FOR TYPE-II AVALANCHES 



A systematic technique for the determination of the avalanche exponents a and v is the so called moment analysis [3| . 
We illustrate how this works for the size distribution of the BTW model. The analysis confirms our result v = 0.49 
for the BTW model for type-II avalanches. 

Let Psizc{s; N) be the PDF of avalanche size s in the two-dimensional BTW model with N'^ sites and a threshold 
{x)/3. Let us assume finite-size scaling (FSS). This means that for iV, s ^ 1 we have 

Psizc(s;iV)cxs-''G(s/iV^), (13) 

for some exponent D > 0, where G{r) falls off quickly for r > 1. From Eq. (fT3l ) we observe that 

/•OO 

{s") oc N^^^+'i-''^ r"-" G{r) dr . 

Jl/NO 

The integral tends to a constant as A^ ^ cxi so we have (s'?) ^ ]\[D{i+q-^) _ Thus, if we plot (s'?) versus A^ in a log-log 
plot the slope of a fitted straight line will give us an estimate of the exponent ({q) = D{1 + q — j^) for q = 1,2,3,.... 
FigjS^ shows the computation of moments of s as a function of N for type-II avalanches obtained from simulation 
of the BTW model, and FiglSjD shows the exponent ({q) versus q. We find that ({q) ~ q"^'^^, hence that D = 2.81. 
Moreover, when plotted in a log-log plot, the intersection of C{q) with the first axis corresponds to the value i/ — 1. 
Hence, based on our predictions this intersection should be in the point 0.49. This value is plotted as a dotted 
horizontal line in FiglSb, confirming our result with good accuracy. 

Strictly, this method requires that we have data collapse when re-scaling the avalanche sizes by s/N^ . The shape 
of the scaling function G{r) can be seen by plotting s'^Psizcis',N) versus s/N^. This is shown in Fig. [6l The data 
collapse for the available system sizes is not perfect, but it is better than for type-I avalanches. In fact, it seems 
that we might have convergence to a Af-independent scaling function as A^ ^- oo, and that the lack of data collapse 
observed here is simply due to the fact that we are not able to simulate sufficiently large sandpiles. Thus the general 
lack of scaling for for type-I avalanches, including the finite-size scaling, seems to be restored for type-II avalanches. 
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FIG. 6: Attempted data collapse for the avalanches in the BTW model with respect to a threshold (X)/3. We have plotted 
s"-^ /r VA-z.^{s') ds' as a function of s/L° with D = 2.81 for TV = 64, 256, 1024, 2048. 
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FIG. 7: The PDF of the fluctuations in the toppling activity of the slowly driven BTW sandpile together with the corresponding 
PDF produced from the stochastic model. The dashed line is a fitted stretched exponential. 



THE PDFS OF THE TOPPLING ACTIVITY 



When calculating the PDFs of the toppUng activity signal we have to distinguish between the weakly and strongly 
driven sandpiles. For the weakly driven sandpiles the toppling activity covers a large range and the high-activity tail 
of the PDF decays like a stretched exponential. This property can be reproduced by simulations of Eq. ^ where 
f{X) decays exponentially to zero as X increases. Fig 17] compares the PDF obtained from such a simulation of the 
stochastic differential equation (run with H = 0.37) with the PDF of the weakly driven BTW model. The stochastic 
model is run by initiating new avalanches (realizations) whenever the previous avalanche terminates, thus representing 
the classical slow drive of the sandpile model. Since the Hurst exponents of the Zhang and BTW models are different 
from 1/2, the application of the Fokker-Planck equation can not be used to calculate the shape of the PDFs, and thus 
we have to rely on numerical simulations. 

It is also interesting to study the PDFs of toppling activity in strongly driven sandpiles. In particular in the light 
of the recent claims that this toppling activity belongs to a class of fluctuating global quantities with a universal 
non-gaussian shape [lOl . Illl | . These PDFs are reminiscent of distributions derived from the extreme value theory of 
statistics |I2l |l3|, which deals with a sequence of n identically distributed, independent random variables yi, . . . , y„. 
If the tail of the PDF of these variables decays faster than any power-law, the PDF of the m'th largest value drawn 
from each of M realizations of this sequence converge to the Gumbel class of stable distributions in the limit M — > oo: 



Gkiy)^K{e 



-?fe-«)e-e""""° 



Y' 



(14) 



The distribution for the largest value in each realization {m = 1) is often called the Fischer-Tippet-Gumbel distribution 
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FIG. 8: The PDF of the fluctuations in the toppling activity of the strongly driven Zhang sandpile driven BTW sandpile for 
different values of A (see text) . The dashed line is a fitted Gaussian. 
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FIG. 9: a) Part of a time series for the toppling activity in a strongly driven Zhang sandpile. b) The increments Ax{t) = 
x{t + At) — x{t) of the signal in (a). As for the slowly driven case, c) The conditional variance of the increments Ax given the 
value of X. d) The conditional mean of increments. This curve is fitted to a parabola. 



(FTG). The FTG with zero mean and unit variance requires K = (_ = tt/^/Q and s — ^VG/tt, where 7 « 0.58 is the 
so-called Euler constant. 

A distribution obtained from the spin-wave approximation to the 2D XY model for equilibrium crititical fluctua- 
tions in a finite-size magnetized system is the so-called Bramwell-Holdsworth-Pinton (BHP) distribution [ll|, which 
corresponds to Eq. p^ with k having the non-integer value m = 7r/2 ss 1.57. With K — 2.16, a — 1.58, £, — 0.93 and 
s = 0.37 this distribution has zero mean and unit variance. The difference between the normalized FTG and BHP 
distributions is rather small, and it is difficuh to distinguish between the two based in experimental and numerical 
data. 

Analysis of our simulations show the claim that the toppling activity in strongly driven sandpiles has a PDF similar 
to the BHP or FTG distributions is wrong, unless the driving rate of the system is fine tuned to some particular 
value. Actually, the normalized PDFs of the toppling activity is only insensitive to variation of the driving rate in 
the limits of weak and strong drive. In the limit of weak drive the PDFs are stretched exponentials as shown in Fig[7| 
and in the limit of strong drive the PDFs are close to Gaussian. 

FigOshows the PDFs of the toppling activity in the strongly driven Zhang model where the feeding rate is a Poisson 
process with characteristic time scales A = (feeding one unit e of mass in every time step), A = 2.3, and A = 3.0. For 
A = and A — 2.3 the sandpile is running, in the meaning that avalanches never terminate. For A = 3, the sandpile 
is no longer running, and we see a deviation from the Gaussian shape in the left tail of the PDF. 

The Gaussian PDFs for strongly driven sandpiles can actually be explained in terms of the stochastic model in 
Eq. ([2|). The idea is that the range of fluctuations is conflned by the drift term f{X), which now has a different shape 
than in the slowly driven sandpile. Fig[5] shows the shape of both the diffusion coefficient D{X) = a'^X/2 and the 
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FIG. 10: a) Part of a time series for the kinetic energy in the 2D turbulence simulation, b) The increments AX{t) = 
X{t + At) — X{t) of the signal in (a). As for the toppling activity in sandpiles we see that the increments are large when the 
kinetic energy itself is large, c) The conditional variance of the increments AX given the value of X. The inset shows the 
logarithm of this variance versus X. The dashed curves corresponds to a fitted exponential function, d) The conditional mean 
of increments. We see that the mean is approximately constant for a large range of X. 



drift term for the strongly driven Zhang model. The diffusion coefficient has the same form as for the slowly driven 
sandpile, whereas the drift term is well approximated by a parabola: f{X) = —aX^ + bX + c. As explained in Sec. U 
this drift term confines the fluctuations in toppling activity to a bounded region around the positive root Xc of /(AT). 
Due to the higher rate of random feeding, the memory effects described by the Hurst exponent H ^ 1/2 becomes 
inessential in the strongly driven sandpile, allowing us to give an approximate description of the time dependent PDF 
through the Fokker-Planck equation 



dP _ d 
~dt '^dX 

Stationary solutions of this equation must satisfy 



f{X)P 



92 



2 dX^ 



(XP). 



(15) 



-^(^W^ 



-—2i^P)=^' 



2 dX 
and substituting —aX^ + bX + c for /(AT), we have the Gaussian solution 



PiX) 



1 



27rS 



exp 



2S2 



with 



^i — — and Zj = —= — . 
a ^/2 a 

To further substantiate that the toppling dynamics of sandpiles is fundamentally different from the fluctuating 
quantities that give rise to BHP and FTG distributions we apply the above analysis to the fluctuations in total 
kinetic energy in a simulation of two-dimensional Navier- Stokes turbulence. The two-dimensional geometry is chosen 
because of the inverse energy cascade in fc-space caused by merging of small vortices and formation of large structures, 
reminiscent of formation of large avalanches from emerging from localized random perturbation in sandpile models. 
In the simulation energy is injected via a source term on a characteristic, small spatial scale throughout the simulation 
area, and is dissipated as loss through the open boundary. 

We compute the conditional mean and conditional variance of the increment process. These results are presented in 
FigHni We observe that contrary to the sandpile models, the diffusion coefficient for this process grows exponentially 
with X . Moreover, the drift term can be approximated by a small positive constant except when the kinetic energy is 
very large. This leads us to the stochastic differential equation Eq. ([3]), and the corresponding Fokker-Planck equation 
is 



dP 
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dP 
"dX 



2 dxA^ 



(16) 
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FIG. 11: The normalized PDF for the fluctuation of kinetic energy in the 2D turbulence simulation together with a normalized 
FTG PDF (smooth solid curve), the normalized BHP PDF (smooth dashed curve) and a PDF obtained from the simulation 
of Eq. [3] (dotted curve). 



:(e2«^p) 



Stationary solutions of this equation must satisfy 

,dP 
dX ' 2 dX^ 

We can put c — 1 without loss of generahty, and obtain the solution 

1 



= 0. 



P(X) = Le-^x-,)/p^-e-i--^y^ 



where /3 = l/2a and /i = (l/2a) log (l/2a). This is the standard FTG distribution. FigfTTj shows the normalized 
PDF for the fluctuations in total kinetic energy obtained from the fluid simulation along with a normalized FTG 
distribution and the PDF obtained from the simulation of Eq. ^ . 



VI. CONCLUSIONS 



When output from simple models for avalanching systems are compared to observational data of real- world systems 
one has to deal with the problem of establishing a correspondence between the model variables and the observables 
of the natural system. In general, this is usually not an obvious task, since the model is usually not derived from first 
physical principles. The observation may be spatiotemporal, or just temporal. Likewise we may choose to analyze 
the spatiotemporal output from an avalanche model, or just some spatially integrated quantity like the total activity 
variable in a sandpile, presented as a time series. In this paper we have focused on the latter, and in particular on 
reproducing the statistical properties of such time-series by modeling the stochastic process by means of a stochastic 
differential equation. 

Since our interest is avalanching dynamics we focus on avalanche statistics, and we have to face the problem of 
how to define what an avalanche is when our observational data is in the form of a time series. In general we cannot 
expect that such a definition is necessarily equivalent with a definition based on spatiotemporal data, at least not in 
those cases where feeding of new "sand grains" occurs while avalanches are running. In those cases several spatially 
separated avalanches may run simultaneously, but these cannot be separated in a pure temporal analysis. For such 
a continuously driven model system the temporal signal may never be zero, and in an observational time series noise 
also makes it impossible to use a zero signal condition to separate an avalanching state from a quiet state. Thus, the 
natural way out is to define the avalanching state by means of a threshold level on the activity signal, giving rise to 
the concept of type-II avalanches. 

In this paper we have shown that the toppling activity in sandpiles, and also global kinetic energy in a 2D fiuid 
simulation, can be modeled by stochastic differential equations. The modeling clarifies that the main discrepancy 
between the mean-field approach and the actual BTW and Zhang models is related to the Hurst exponent of the 
activity process, which is different for the two sandpile models. It also clarifies the origin of the differences between the 
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scaling exponents for type-I and type-II avalanches, and why type-II avalanches exhibit clearer scaling characteristics 
than type-I avalanches. 

It follows from the theory how to rescale coordinates to approach the thermodynamic limit, and the results obtained 
for finite-size scaling in the BTW model in Sec. II VI gives a strong indication that this limit actually exists. 

For continuously driven sandpiles the stochastic equation can be cast into a Fokker-Planck equation due to the loss 
of memory caused by the random feeding. This allows an analytic solution for the activity PDF, which is a Gaussian 
for the sandpile models, but gives rise to the FTG distribution for the 2D fluid simulation. These results show that 
the non-Gaussian universal PDF described in |10l . Ill| is not relevant for strongly driven sandpiles, but may be so 
for certain turbulence models. The stochastic theory relates the difference between Gaussian and FTG distributed 
activity signal to a the difference between a linear and exponential X dependence of the diffusion coefficient in the 
stochastic equation. 
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